	subroutine rinj(ireg)
        USE vast_kind_param
        use corgan_com_M 
        use cindex_com_M 
        use numpar_com_M 
        use cprplt_com_M
        use cophys_com_M 
        use ctemp_com_M
        use blcom_com_M
        use objects_com_M  
	implicit none
	logical :: verbose
	integer :: ip,ncellsp2,n,npinj,np,npcel,ireg
        integer,dimension(itdim) :: ijkctmp2
	real(8) :: sbottom,sutop,delz,vtherm,flxrnd,flxdir,fluxin,  &
         qbottom,qtop,qout,qinj,diepar,vmag,ws,pxi2,peta2,pzta2,&
         th,wi,wip,wj,wjp,wk,wkp,uinj,vinj,winj,deninj, delx,shaperx,shapery,segno,derf
      real(double) , external :: ranf,generate
!
!	A subroutine to inject particles
!
	verbose=.false.
        pi=acos(-1.)
        npcel=npcelx(ireg)*npcely(ireg)*npcelz(ireg)
	segno = sign(1.,qom(icoi(ireg)))

        if(inject_end(ireg)) then
!
!	END Minimum Z boundary
!
        call celindex(2,ibp1,2,jbp1,2,2,iwid,jwid,kwid,ijkctmp,  &
         ijkctmp,ncellsp,ijkctmp2,ncellsp2,ncellsp2)

        do n=1,ncellsp
            ijk=ijkctmp(n)
              call boundary_profile(ireg,vtherm,uinj,vinj,winj,deninj,1)	
	      flxrnd=exp(-winj**2/vtherm**2)*vtherm/2./sqrt(pi)
	      flxdir=.5*winj*(1.+derf(winj/vtherm))
	      fluxin=deninj*(flxrnd+flxdir)
            delz=z(ijk+kwid)-z(ijk)
            npinj=int(fluxin*dt*npcel/rhr(ireg)/delz)
	    if(npinj.eq.0) then
!	   if(ranf().lt.fluxin*dt*npcel/deninj/delz) then
!	      npinj=1
!	      qinj=deninj*vol(ijk)/npcel
!           endif
	    npinj=1
            qinj=fluxin*vol(ijk)/delz*dt
            else
	    qinj=fluxin*vol(ijk)/delz*dt
	end if

	qinj=qinj*sign(1.d0,qom(ireg))

        if (verbose) write(*,*) 'rinj:  npinj=',npinj
	do ip=1,npinj
	np=iphead(1)
        if (verbose) write(*,*) 'rinj: np,link(np)=',np,link(np)

	if(np.gt.0) then

 	iphead(1)=link(np)
	link(np)=iphd2(ijk)
	iphd2(ijk)=np

        pxi2=ranf()
        peta2=ranf()
        pzta2=0d0
!
!     calculate interpolation weights
!
         wi = 1. - pxi2
         wip = pxi2
!
         wj = 1. - peta2
         wjp = peta2
!
         wk = 1. - pzta2
         wkp = pzta2
!
         wate(ijk,1) = wip*wj*wk
         wate(ijk,2) = wip*wjp*wk
         wate(ijk,3) = wi*wjp*wk
         wate(ijk,4) = wi*wj*wk
!
         wate(ijk,5) = wip*wj*wkp
         wate(ijk,6) = wip*wjp*wkp
         wate(ijk,7) = wi*wjp*wkp
         wate(ijk,8) = wi*wj*wkp


         px(np) = (wate(ijk,1)*x(ijk+iwid)   &
            +(wate(ijk,2)*x(ijk+iwid+jwid)&
            +(wate(ijk,3)*x(ijk+jwid)  &
            +(wate(ijk,4)*x(ijk)  &
            +(wate(ijk,5)*x(ijk+iwid+kwid)  &
            +(wate(ijk,6)*x(ijk+iwid+jwid+kwid)   &  
            +(wate(ijk,7)*x(ijk+jwid+kwid)&
            +(wate(ijk,8)*x(ijk+kwid)))))))))
         py(np) = (wate(ijk,1)*y(ijk+iwid)   &
            +(wate(ijk,2)*y(ijk+iwid+jwid)&
            +(wate(ijk,3)*y(ijk+jwid)  &
            +(wate(ijk,4)*y(ijk)  &
            +(wate(ijk,5)*y(ijk+iwid+kwid)  &
            +(wate(ijk,6)*y(ijk+iwid+jwid+kwid)   &  
            +(wate(ijk,7)*y(ijk+jwid+kwid)&
            +(wate(ijk,8)*y(ijk+kwid)))))))))
         pz(np) = (wate(ijk,1)*z(ijk+iwid)   &
            +(wate(ijk,2)*z(ijk+iwid+jwid)&
            +(wate(ijk,3)*z(ijk+jwid)  &
            +(wate(ijk,4)*z(ijk)  &
            +(wate(ijk,5)*z(ijk+iwid+kwid)  &
            +(wate(ijk,6)*z(ijk+iwid+jwid+kwid)   &  
            +(wate(ijk,7)*z(ijk+jwid+kwid)&
            +(wate(ijk,8)*z(ijk+kwid)))))))))

	qpar(np)=qinj/npinj
	ico(np)=icoi(ireg)
        vmag = sqrt(-2.*log(1.-.999999*ranf()))
        th = 2.*pi*ranf()
	!shaperx = tanh((z(ijk)-0.5*(zf+ze))/el_rec)/cosh((z(ijk)-0.5*(zf+ze))/el_rec) 
	shaperx=1d0
        up(np)=segno*uinj*shaperx+vmag*siepx(ireg)*cos(th)
        vp(np)=segno*vinj+vmag*siepy(ireg)*sin(th)
	wp(np)=generate(winj,siepz(ireg)*sqrt(2.))
	pz(np)=pz(np)+dt*ranf()*wp(np)
!      vmag = sqrt(-2.*log(1.-.999999*ranf()))
!        th = 2.*pi*ranf()
!        wp(np)=segno*winj*0d0+vmag*siepz(ireg)*sin(th)
!	pz(np)=pz(np)+dt*ranf()*winj	
	else

 	write(*,*)'WARNING: Not Enough Particles'

	end if
	
	enddo

	enddo
        endif !inject_end


        if(inject_front(ireg)) then
!
!	Front
!
        call celindex(2,ibp1,2,jbp1,kbp1,kbp1,iwid,jwid,kwid,ijkctmp,  &
         ijkctmp,ncellsp,ijkctmp2,ncellsp2,ncellsp2)

        do n=1,ncellsp
            ijk=ijkctmp(n)
              call boundary_profile(ireg,vtherm,uinj,vinj,winj,deninj,-1)	
	      flxrnd=exp(-winj**2/vtherm**2)*vtherm/2./sqrt(pi)
	      flxdir=.5*winj*(1.+derf(winj/vtherm))
	      fluxin=deninj*(flxrnd+flxdir)
            delz=z(ijk)-z(ijk-kwid)
        npinj=int(fluxin*dt*npcel/rhr(ireg)/delz)
	if(npinj.eq.0) then
!	   if(ranf().lt.fluxin*dt*npcel/deninj/delz) then
!	      npinj=1
!	      qinj=deninj*vol(ijk)/npcel
!           endif
	npinj=1
        qinj=fluxin*vol(ijk)/delz*dt
        else
	    qinj=fluxin*vol(ijk)/delz*dt
	end if

	qinj=qinj*sign(1.d0,qom(ireg))

	if(verbose) write(*,*) 'rinj: Front npinj',npinj
	do ip=1,npinj
	np=iphead(1)

	if(np.gt.0) then

 	iphead(1)=link(np)
	link(np)=iphd2(ijk)
	iphd2(ijk)=np

        pxi2=ranf()
        peta2=ranf()
        pzta2=1.d0

!
!     calculate interpolation weights
!
         wi = 1. - pxi2
         wip = pxi2
!
         wj = 1. - peta2
         wjp = peta2
!
         wk = 1. - pzta2
         wkp = pzta2
!
         wate(ijk,1) = wip*wj*wk
         wate(ijk,2) = wip*wjp*wk
         wate(ijk,3) = wi*wjp*wk
         wate(ijk,4) = wi*wj*wk
!
         wate(ijk,5) = wip*wj*wkp
         wate(ijk,6) = wip*wjp*wkp
         wate(ijk,7) = wi*wjp*wkp
         wate(ijk,8) = wi*wj*wkp

         px(np) = (wate(ijk,1)*x(ijk+iwid)   &
            +(wate(ijk,2)*x(ijk+iwid+jwid)&
            +(wate(ijk,3)*x(ijk+jwid)  &
            +(wate(ijk,4)*x(ijk)  &
            +(wate(ijk,5)*x(ijk+iwid+kwid)  &
            +(wate(ijk,6)*x(ijk+iwid+jwid+kwid)   &  
            +(wate(ijk,7)*x(ijk+jwid+kwid)&
            +(wate(ijk,8)*x(ijk+kwid)))))))))
         py(np) = (wate(ijk,1)*y(ijk+iwid)   &
            +(wate(ijk,2)*y(ijk+iwid+jwid)&
            +(wate(ijk,3)*y(ijk+jwid)  &
            +(wate(ijk,4)*y(ijk)  &
            +(wate(ijk,5)*y(ijk+iwid+kwid)  &
            +(wate(ijk,6)*y(ijk+iwid+jwid+kwid)   &  
            +(wate(ijk,7)*y(ijk+jwid+kwid)&
            +(wate(ijk,8)*y(ijk+kwid)))))))))
         pz(np) = (wate(ijk,1)*z(ijk+iwid)   &
            +(wate(ijk,2)*z(ijk+iwid+jwid)&
            +(wate(ijk,3)*z(ijk+jwid)  &
            +(wate(ijk,4)*z(ijk)  &
            +(wate(ijk,5)*z(ijk+iwid+kwid)  &
            +(wate(ijk,6)*z(ijk+iwid+jwid+kwid)   &  
            +(wate(ijk,7)*z(ijk+jwid+kwid)&
            +(wate(ijk,8)*z(ijk+kwid)))))))))

	qpar(np)=qinj/npinj
	ico(np)=icoi(ireg)
        vmag = sqrt(-2.*log(1.-.999999*ranf()))
        th = 2.*pi*ranf()
	!shaperx = tanh((z(ijk)-0.5*(zf+ze))/el_rec)/cosh((z(ijk)-0.5*(zf+ze))/el_rec) 
	shaperx=1d0

       up(np)=segno*uinj*shaperx+vmag*siepx(ireg)*cos(th)
        vp(np)=segno*vinj+vmag*siepy(ireg)*sin(th)
	wp(np)=-generate(winj,siepz(ireg)*sqrt(2.))
	pz(np)=pz(np)+dt*ranf()*wp(np)
!      vmag = sqrt(-2.*log(1.-.999999*ranf()))
!        th = 2.*pi*ranf()
!        wp(np)=segno*winj*0d0+vmag*siepz(ireg)*sin(th)
!	pz(np)=pz(np)+dt*ranf()*winj	
	else

 	write(*,*)'WARNING: Not Enough Particles'

	end if
	
	enddo

	enddo

        endif  !  inject_front

        if(inject_left(ireg)) then
!
!	LEFT
!
        call celindex(2,2,2,jbp1,2,kbp1,iwid,jwid,kwid,ijkctmp,  &
         ijkctmp,ncellsp,ijkctmp2,ncellsp2,ncellsp2)

        do n=1,ncellsp
            ijk=ijkctmp(n)
              call maxwellian_inj(ireg,vtherm,uinj,vinj,winj,deninj,1)	
	      flxrnd=exp(-uinj**2/vtherm**2)*vtherm/2./sqrt(pi)
	      flxdir=.5*uinj*(1.+derf(uinj/vtherm))
	      fluxin=deninj*(flxrnd+flxdir)
            delx=x(ijk+iwid)-x(ijk)
            npinj=int(fluxin*dt*npcel/rhr(ireg)/delx)
	    if(npinj.eq.0) then
!	   if(ranf().lt.fluxin*dt*npcel/deninj/delx) then
!	      npinj=1
!	      qinj=deninj*vol(ijk)/npcel
!           endif
	    npinj=1
            qinj=fluxin*vol(ijk)/delx*dt
            else
	    qinj=fluxin*vol(ijk)/delx*dt
	end if

	qinj=qinj*sign(1.d0,qom(ireg))

        if (verbose) write(*,*) 'rinj:  npinj=',npinj
	do ip=1,npinj
	np=iphead(1)
        if (verbose) write(*,*) 'rinj: np,link(np)=',np,link(np)

	if(np.gt.0) then

 	iphead(1)=link(np)
	link(np)=iphd2(ijk)
	iphd2(ijk)=np

        pxi2=0d0
        peta2=ranf()
        pzta2=ranf() 
!
!     calculate interpolation weights
!
         wi = 1. - pxi2
         wip = pxi2
!
         wj = 1. - peta2
         wjp = peta2
!
         wk = 1. - pzta2
         wkp = pzta2
!
         wate(ijk,1) = wip*wj*wk
         wate(ijk,2) = wip*wjp*wk
         wate(ijk,3) = wi*wjp*wk
         wate(ijk,4) = wi*wj*wk
!
         wate(ijk,5) = wip*wj*wkp
         wate(ijk,6) = wip*wjp*wkp
         wate(ijk,7) = wi*wjp*wkp
         wate(ijk,8) = wi*wj*wkp


         px(np) = (wate(ijk,1)*x(ijk+iwid)   &
            +(wate(ijk,2)*x(ijk+iwid+jwid)&
            +(wate(ijk,3)*x(ijk+jwid)  &
            +(wate(ijk,4)*x(ijk)  &
            +(wate(ijk,5)*x(ijk+iwid+kwid)  &
            +(wate(ijk,6)*x(ijk+iwid+jwid+kwid)   &  
            +(wate(ijk,7)*x(ijk+jwid+kwid)&
            +(wate(ijk,8)*x(ijk+kwid)))))))))
         py(np) = (wate(ijk,1)*y(ijk+iwid)   &
            +(wate(ijk,2)*y(ijk+iwid+jwid)&
            +(wate(ijk,3)*y(ijk+jwid)  &
            +(wate(ijk,4)*y(ijk)  &
            +(wate(ijk,5)*y(ijk+iwid+kwid)  &
            +(wate(ijk,6)*y(ijk+iwid+jwid+kwid)   &  
            +(wate(ijk,7)*y(ijk+jwid+kwid)&
            +(wate(ijk,8)*y(ijk+kwid)))))))))
         pz(np) = (wate(ijk,1)*z(ijk+iwid)   &
            +(wate(ijk,2)*z(ijk+iwid+jwid)&
            +(wate(ijk,3)*z(ijk+jwid)  &
            +(wate(ijk,4)*z(ijk)  &
            +(wate(ijk,5)*z(ijk+iwid+kwid)  &
            +(wate(ijk,6)*z(ijk+iwid+jwid+kwid)   &  
            +(wate(ijk,7)*z(ijk+jwid+kwid)&
            +(wate(ijk,8)*z(ijk+kwid)))))))))

	qpar(np)=qinj/npinj
	ico(np)=icoi(ireg)
        vmag = sqrt(-2.*log(1.-.999999*ranf()))
        ws=vmag*siepz(ireg)
        th = 2.*pi*ranf()
shapery = 1.d0/cosh((pz(np)-0.5*(zf+ze))/el_rec)**2
shapery=1d0
        wp(np)=(winj*segno+ws*cos(th))
        vp(np)=(segno*vinj*shapery+ws*sin(th))
	up(np)=generate(uinj,siepx(ireg)*sqrt(2.))
	px(np)=px(np)+ranf()*dt*up(np)

	else

 	write(*,*)'WARNING: rhrNot Enough Particles'

	end if
	
	enddo

	enddo
        endif !inject_left


if(inject_right(ireg)) then
!
!	RIGHT
!segno
        call celindex(ibp1,ibp1,2,jbp1,2,kbp1,iwid,jwid,kwid,ijkctmp,  &
         ijkctmp,ncellsp,ijkctmp2,ncellsp2,ncellsp2)

        do n=1,ncellsp
            ijk=ijkctmp(n)
              call maxwellian_inj(ireg,vtherm,uinj,vinj,winj,deninj,-1)	
	      flxrnd=exp(-uinj**2/vtherm**2)*vtherm/2./sqrt(pi)
	      flxdir=.5*uinj*(1.+derf(uinj/vtherm))
	      fluxin=deninj*(flxrnd+flxdir)
            delx=x(ijk+iwid)-x(ijk)
            npinj=int(fluxin*dt*npcel/rhr(ireg)/delx)
	    if(npinj.eq.0) then
!	   if(ranf().lt.fluxin*dt*npcel/deninj/delx) then
!	      npinj=1
!	      qinj=deninj*vol(ijk)/npcel
!           endif
	    npinj=1
            qinj=fluxin*vol(ijk)/delx*dt
            else
	    qinj=fluxin*vol(ijk)/delx*dt
	end if

	qinj=qinj*sign(1.d0,qom(ireg))

        if (verbose) write(*,*) 'rinj:  npinj=',npinj
	do ip=1,npinj
	np=iphead(1)
        if (verbose) write(*,*) 'rinj: np,link(np)=',np,link(np)

	if(np.gt.0) then

 	iphead(1)=link(np)
	link(np)=iphd2(ijk)
	iphd2(ijk)=np

        pxi2=1d0
        peta2=ranf()
        pzta2=ranf()
!
!     calculate interpolation weights
!
         wi = 1. - pxi2
         wip = pxi2
!
         wj = 1. - peta2
         wjp = peta2
!
         wk = 1. - pzta2
         wkp = pzta2
!
         wate(ijk,1) = wip*wj*wk
         wate(ijk,2) = wip*wjp*wk
         wate(ijk,3) = wi*wjp*wk
         wate(ijk,4) = wi*wj*wk
!
         wate(ijk,5) = wip*wj*wkp
         wate(ijk,6) = wip*wjp*wkp
         wate(ijk,7) = wi*wjp*wkp
         wate(ijk,8) = wi*wj*wkp


         px(np) = (wate(ijk,1)*x(ijk+iwid)   &
            +(wate(ijk,2)*x(ijk+iwid+jwid)&
            +(wate(ijk,3)*x(ijk+jwid)  &
            +(wate(ijk,4)*x(ijk)  &
            +(wate(ijk,5)*x(ijk+iwid+kwid)  &
            +(wate(ijk,6)*x(ijk+iwid+jwid+kwid)   &  
            +(wate(ijk,7)*x(ijk+jwid+kwid)&
            +(wate(ijk,8)*x(ijk+kwid)))))))))
         py(np) = (wate(ijk,1)*y(ijk+iwid)   &
            +(wate(ijk,2)*y(ijk+iwid+jwid)&
            +(wate(ijk,3)*y(ijk+jwid)  &
            +(wate(ijk,4)*y(ijk)  &
            +(wate(ijk,5)*y(ijk+iwid+kwid)  &
            +(wate(ijk,6)*y(ijk+iwid+jwid+kwid)   &  
            +(wate(ijk,7)*y(ijk+jwid+kwid)&
            +(wate(ijk,8)*y(ijk+kwid)))))))))
         pz(np) = (wate(ijk,1)*z(ijk+iwid)   &
            +(wate(ijk,2)*z(ijk+iwid+jwid)&
            +(wate(ijk,3)*z(ijk+jwid)  &
            +(wate(ijk,4)*z(ijk)  &
            +(wate(ijk,5)*z(ijk+iwid+kwid)  &
            +(wate(ijk,6)*z(ijk+iwid+jwid+kwid)   &  
            +(wate(ijk,7)*z(ijk+jwid+kwid)&
            +(wate(ijk,8)*z(ijk+kwid)))))))))

	qpar(np)=qinj/npinj
	ico(np)=icoi(ireg)
        vmag = sqrt(-2.*log(1.-.999999*ranf()))
        ws=vmag*siepz(ireg)
        th = 2.*pi*ranf()
shapery = 1.d0/cosh((pz(np)-0.5*(zf+ze))/el_rec)**2
shapery=1d0
        wp(np)=(winj*segno+ws*cos(th))
        vp(np)=(segno*vinj*shapery+ws*sin(th))
	up(np)=-generate(uinj,siepx(ireg)*sqrt(2.))
	px(np)=px(np)+ranf()*dt*up(np)

	else

 	write(*,*)'WARNING: Not Enough Particles'

	end if
	
	enddo

	enddo
        endif !inject_right

	return
	end

	subroutine boundary_profile(ireg,vtherm,uinj,vinj,winj,deninj,flag)
        USE vast_kind_param
        use corgan_com_M 
        use cindex_com_M 
        use numpar_com_M 
        use cprplt_com_M
        use cophys_com_M 
        use ctemp_com_M
        use blcom_com_M
        use objects_com_M  
	implicit none
	integer :: flag,ireg
	real(8) :: a_m,o_m,g_x,v_int,bavg2,vavg,timescale,z_c
	real(8) :: vtherm,uinj,vinj,winj,deninj,zone,ztwo,fx,harris,segno
!	flag is scalar product of unit verctor with inward normal 
!	flag=1  lower end
!	flag=-1 upper end

                     ipjk = ijk + iwid 
                     ipjpk = ijk + iwid + jwid 
                     ijpk = ijk + jwid 
!
                     ijkp = ijk + kwid 
                     ipjkp = ijk + iwid + kwid 
                     ijpkp = ijk + jwid + kwid 
                     ipjpkp = ijk + iwid + jwid + kwid 
!
                     zone = 0.25*(z(ipjk)+z(ipjpk)+z(ijpk)+z(ijk)) 
                     ztwo = 0.25*(z(ipjkp)+z(ipjpkp)+z(ijpkp)+z(ijkp)) 
!
!
!       BIRN profile
!
                        fx = (1. + eps_rec/rnu_rec*x(ijk)/el_rec)**(-rnu_rec) 
!
!       Quasiparabolic profile (Pritchett)
!
!        fx=exp(-x(ijk)*eps_rec/el_rec)
!
    harris = el_rec*fx/(ztwo - zone)*(tanh((ztwo - 0.5*(zf&
             + ze))*fx/el_rec) - tanh((zone - 0.5*(zf + ze))*fx/el_rec)) 

	harris=harris+rhr(nrg)
	
	timescale=bx_ini
	a_m=ey_ext
	o_m=.05d0*timescale

	z_c=kbar*dz/2.d0
	g_x=sin(pi*x(ijk)/dx/dble(ibar))**2
	v_int=2.d0*a_m*o_m*tanh(o_m*t)/cosh(o_m*t)**2
	if(shock) then
	    g_x=1d0
        v_int=ey_ext
    end if
	winj=v_int*g_x
!    segno = sign(1.,qom(icoi(ireg)))
	!write(*,*)'INFLOW SPEED=',vavg

    vtherm=siepz(ireg)*sqrt(2.)

!    if(abs(ey_ext).lt.1d-10) then
!        winj=wvi(ireg)*segno
!    else
!        winj=vavg
!	end if 
!    call vector_product_norm(exavp(ijk+flag*iwid),eyavp(ijk+flag*iwid),ezavp(ijk+flag*iwid), &
!		                 bxv(ijk+flag*iwid),byv(ijk+flag*iwid),bzv(ijk+flag*iwid), &
!						 uinj,vinj,winj)
	deninj=harris
!    deninj=rhr(ireg)
    uinj=uinj*0d0+uvi(ireg)
    vinj=vinj*0d0+vvi(ireg)
    winj=winj+wvi(ireg)
	
!	write(*,*)'injkector',ireg,winj,deninj

    end subroutine boundary_profile

	subroutine maxwellian_inj(ireg,vtherm,uinj,vinj,winj,deninj,flag)
        USE vast_kind_param
        use corgan_com_M 
        use cindex_com_M 
        use numpar_com_M 
        use cprplt_com_M
        use cophys_com_M 	
        use ctemp_com_M
        use blcom_com_M
        use objects_com_M  
	implicit none
	integer :: flag,ireg
	real(double) :: a_m,o_m,g_x,v_int,bavg2,vavg,timescale,z_c
	real(double) :: vtherm,uinj,vinj,winj,deninj,zone,ztwo,fx,shaperx,segno, harris,xc,zc
	real(double) , external :: ranf
!	flag is scalar product of unit verctor with inward normal 
!	flag=1  lower end
!	flag=-1 upper end
!
                     ipjk = ijk + iwid 
                     ipjpk = ijk + iwid + jwid 
                     ijpk = ijk + jwid 
!
                     ijkp = ijk + kwid 
                     ipjkp = ijk + iwid + kwid 
                     ijpkp = ijk + jwid + kwid 
                     ipjpkp = ijk + iwid + jwid + kwid 
!
                     xc = 0.125*(x(ipjk)+x(ipjpk)+x(ijpk)+x(ijk)+x(ipjkp)+x(ipjpkp)+x&
                          (ijpkp)+x(ijkp)) 
                     zc = 0.125*(z(ipjk)+z(ipjpk)+z(ijpk)+z(ijk)+z(ipjkp)+z(ipjpkp)+z&
                          (ijpkp)+z(ijkp)) 
!

    segno = sign(1.,qom(icoi(ireg)))
    

    if (forcefree) then
        harris=1d0
		shaperx = tanh((zc-0.5*(zf+ze))/el_rec)/cosh((zc-0.5*(zf+ze))/el_rec) 
    else
        harris = (1.-epsil**2)/ &
                 (epsil*cos((xc-.5*(xr+xl))/el_rec) + cosh((zc-.5*(zf+ze))/el_rec))**2
! density perturbation that is consistent with initial insland perturbation
!        harris = harris - pert_gem*bx_ini*((2.d0*pi/(xr-xl))**2 + &
!             (pi/(zf-ze))**2)*cos(2.d0*pi*(xc-0.5d0*(xr+xl))/(xr-xl)) &
!             *cos(pi*(zc-.5d0*(zf+ze))/(zf-ze))        
        shaperx=1d0
	end if
    vtherm=siepx(ireg)*sqrt(2.)
    !uinj=uvi(ireg)*segno*shaperx
    call vector_product_norm(exavp(ijk+flag*iwid),eyavp(ijk+flag*iwid),ezavp(ijk+flag*iwid), &
		                 bxv(ijk+flag*iwid),byv(ijk+flag*iwid),bzv(ijk+flag*iwid), &
						 uinj,vinj,winj)
    deninj=harris+rhr(nrg)
		  
    !deninj=abs(qdnv(ijk+flag*iwid,ireg))

    uinj=flag*(uinj+uvi(ireg))
    vinj=vinj*0d0+vvi(ireg)
    winj=winj*0d0+wvi(ireg)

	end subroutine maxwellian_inj
	
